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ABSTRACT 

We present evidence that low-mass starless cores, the simplest units of star formation, are systemat¬ 
ically differentiated in their chemical composition. Molecules including CO and CS almost vanish near 
the core centers, where the abundance decreases by at least one or two orders of magnitude with respect 
to the value in the outer core. At the same time, the N 2 H“'' molecule has a constant abundance, and the 
fraction of NH 3 increases toward the core center. Our conclusions are based on a systematic study of 
5 mostly-round starless cores (L1498, LI495, L1400K, L1517B, and LI544), which we have mapped in 
C^®O(I-0), CS(2-1), N 2 H+(I- 0 ), NH 3 ( 1 , 1 ) and (2,2) and the 1.2 mm continuum (complemented with 
C^^O(l-O) and C^'^S(2-1) data for some systems). For each core we have built a spherically symmet¬ 
ric model in which the density is derived from the 1.2 mm continuum, the kinetic temperature from 
NH 3 , and the abundance of each molecule is derived using a Monte Carlo radiative transfer code which 
simultaneously fits the shape of the central spectrum and the radial profile of integrated intensity. Re¬ 
garding the cores for which we have C^^O(l-O) and C^‘*S(2-1) data, the model fits these observations 
automatically when the standard isotopomer ratio is assumed. 

As a result of this modeling, we also find that the gas kinetic temperature in each core is constant 
at approximately 10 K. In agreement with previous work, we find that if the dust temperature is also 
constant, then the density profiles are centrally flattened, and we can model them with a single analytic 
expression. We also find that for each core the turbulent linewidth seems constant in the inner 0.1 pc. 

The very strong abundance drop of CO and CS toward the center of each core is naturally explained by 
the depletion of these molecules onto dust grains at densities of 2-6 xlO^ cm“^. N 2 H''' seems unaffected 
by this process up to densities of several 10® or even 10® cm“®, while the NH 3 abundance may be 
enhanced by its lack of depletion and reactions triggered by the disappearance of CO from the gas 
phase. 

With the help of the Monte Carlo modeling, we show that chemical differentiation automatically ex¬ 
plains the discrepancy between the sizes of CS and NH 3 maps, a problem which has remained unexplained 
for more than a decade. Our models, in addition, show that a combination of radiative transfer effects 
can give rise to the previously observed discrepancy in the linewidth of these two tracers. Although this 
discrepancy has been traditionally interpreted as resulting from a systematic increase of the turbulent 
linewidth with radius, our models show that it can arise in conditions of constant gas turbulence. 

Subject headings: ISM: abundances— ISM: clouds—ISM: molecules—stars: formation 


1. INTRODUCTION 

Dense molecular cores are the basic units of star for¬ 
mation in nearby clouds like Taurus and Perseus, where 
stars like our Sun have been forming over the last few 
million years (e.g., Myers 1999). The study of the phys¬ 
ical structure and kinematics of these cores is therefore 
crucial for our understanding of the star formation pro¬ 
cess, and molecular lines play a role in every step of this 


work. They probe density and temperature through their 
excitation, and turbulent and systematic motions through 
their linewidth and Doppler shifts. For this reason, chem¬ 
ical anomalies in the core gas can hinder our attempt to 
understand core properties, as the lack of a full theory 
of core chemical composition has made it a standard prac¬ 
tice to assume a homogeneous abundance for all molecular 
species. 

The presence of chemical inhomogeneities in the star- 
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forming material at scales of dark clouds has been known 
for some time, with TMC-1 and L134N being the most 
studied examples (e.g., Little et al. 1979; Pratap et al. 
1997; Swade 1989; Dickens et al. 2000). The large-scale 
abundance gradients in these clouds seem best explained 
with time-dependent, gas-phase chemistry, implying that 
different condensations have evolved with different time 
scales (see Langer et al. 2000 and van Dishoeck & Blake 
1998 for recent reviews). At the smaller size of the dense 
cores, a series of recent observations has shown that in 
some cases, the abundance of molecules like CO and CS de¬ 
creases toward the core center (L1498: Kuiper et al. 1996; 
Willacy et al. 1998, IC5146: Kramer et al. 1999; Bergin 
et al. 2001, L977: Alves et al. 1999, L1544: Caselli 
et al. 1999, L1689B: Jessop & Ward-Thompson 2001). 
These decreases in abundance have been interpreted as re¬ 
sulting from the depletion of molecules onto dust grains 
at the high densities and cold temperatures occurring in 
dense core interiors (e.g., Bergin & Langer 1997; Charn- 
ley 1997). It is not clear, however, whether these drops in 
abundance are typical of all dense cores, or they are lim¬ 
ited to a small number of objects. To answer this question, 
we have carried out a systematic study of a sample of 5 
starless cores by observing them in a similar manner and 
analyzing their emission with the same radiative transfer 
modeling. 

2. OBSERVATIONS 

Our core sample is listed in Table 1 and was selected 
by inspecting the N 2 H+(J=l- 0 ) maps of Caselli et al. 
( 2001 a) and choosing cores with bright emission and nearly 
round contours. Our goal was to select cores that had the 
simplest internal structure and lacked obvious outside per¬ 
turbation. In order to achieve high spatial resolution, the 
cores were chosen from the nearby Taurus complex (esti¬ 
mated distance of 140 pc, Elias 1978), with the exception 
of L1400K, which is at an estimated distance of 170 pc 
(Snell 1981). Later mapping showed that L1400K devi¬ 
ates significantly from spherical symmetry, which together 
with its non-Taurus membership made a strong case for 
eliminating it from the sample. We decided however to 
retain this source in our sample to avoid any possible bias. 

As we will show below, this core presents the same chem¬ 
ical behavior as the others, but its elongated shape makes 
it the hardest and most uncertain to model. 

The starless classification of the cores in our sample is 
based on the analysis by Benson & Myers (1989), who 
found that they are not associated with IRAS or NIR 
sources. Unlike Benson & Myers (1989), however, we 
consider L1495 to be starless, given the large separation 
between this core and IRAS 04112-1-2803 (aka CW Tau), 
and the lack of apparent connection between the two. Ac¬ 
cording to Myers et al. (1987), the lack of IRAS detection 
in a Taurus core implies a luminosity limit of < 0.1 Lq for 
a possible embedded source. 

We observed our sample of cores in NH 3 with the 100m 
telescope of the MPIfR in Effelsberg in 1998 October and 
2001 May. The simultaneous (J,K)=(1,1) and (2,2) ob¬ 
servations were done in frequency switching (FSW) mode 
with a throw of 4 MHz, and used the AK90 autocorrela- 

^ FCRAO is supported in part by the National Science Foundation 
Metropolitan District Commission, Commonwealth of Massachusetts. 


tor to provide a velocity resolution of 0.03 or 0.06 km s“^, 
depending on configuration. Cross scans on continuum 
point sources were used to estimate and correct pointing 
errors, and line observations of L1551 and S140 were used 
to calibrate the data, using as a reference the intensities 
reported by Menten & Walmsley (1985) and Ungerechts 
et al. (1986). The telescope beam size at the observing 
frequencies is approximately 40" (FWHM). 

We observed the same 5 cores in C^®O(J=l-0), 
CS(J=2-1), and N 2 H+(J=l- 0 ) with the FCRAO^ 14m 
telescope in 1999 November, 2000 April, and 2001 April. 
To obtain optically thin tracers for testing our radia¬ 
tive transfer models, L1498, L1517B, and L1544 were 
also observed in C^’^O(J=l-0), and L1544 was observed 
C^‘*S(J=2-1). All observations were done with the SE¬ 
QUOIA array receiver in FSW mode, using frequency 
throws of 8 MHz for N 2 H+(l- 0 ), and 4 MHz for the rest 
of the lines. This resulted in a velocity resolution of 0.06 
km s“^ for N 2 H+(l- 0 ), and 0.03 km s“^ for the rest of 
the lines. The pointing was checked and corrected using 
observations of SiO masers, and the data were converted 
into main beam temperature scale using an efficiency of 
0.55 (Ladd & Heyer 1996). The telescope beam size at 
the frequencies of observation is approximately 50". 

Finally, we observed L1498, L1495, L1400K, and 
L1517B in the 1.2mm continuum with the IRAM 30m tele¬ 
scope in 1999 December. The observations were done in 
on-the-fly mode with the MPIfR 37-channel bolometer ar¬ 
ray (Kreysa et al. 1998), using a scanning speed of 4" s“^, 
a wobbler period of 0.5 s, and wobbler throws of 41" and 
53". Three maps of L1495 were done and later combined, 
and the other sources were observed making single maps. 
The sky optical depth was estimated from sky dips before 
and after each individual map. All data were reduced us¬ 
ing the NIC software (Broguiere et al. 1996), and a global 
calibration factor of 15000 counts per Jy beam“^ was es¬ 
timated from an observation of Uranus. The bolometer 
central frequency is 240 GHz and the telescope beam size 
is approximately 11 ". 

For the analysis of the very narrow lines presented in 
this paper it is necessary to use an accurate set of rest fre¬ 
quencies. Currently available line catalogs (e.g., Pickett et 
al. 1998) do not always have the necessary precision, so a 
combination of accurate laboratory measurements and as¬ 
tronomical observations is still necessary. Here we use the 
set of frequencies recommended by Lee et al. (2001), who 
combine new (unpublished) laboratory determinations by 
C. Gottlieb with astronomical observations of the narrow¬ 
line core L1512, which is an improvement of the frequency 
set used by Tafalla et al. (1998). Following Lee et 
al. (2001), we use these frequencies: 93176.258 MHz 
for N 2 H+(JFiF= 101 - 012 ), 97980.953 MHz for CS(J=2- 
1), 109782.172 for Ci®O(J=l-0), and 23694.4949 MHz for 
NH3(J,K=1,1). 

3. RESULTS 

Figure 1 shows our full data set in the form of integrated 
intensity maps. For each starless core (L1498, L1495, 
L1400K, L1517B, and L1544), the three top panels dis¬ 
play maps of what are expected to be column density trac- 
under grant AST 94-20159, and is operated with permission of the 
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Fig. 1.— Maps of 1.2mm continuum (IRAM 30m telescope), C^®O(l-0), C^^O(l-O), CS(2-1), N2H+(l-0) (FCRAO telescope) and NHsfl,!) 
(100m telescope) for L1498, L1495, L1400K, L1517B, and L1544 (L1544 1.2mm continuum map from Ward-Thompson et al. 1999). Central 
coordinates are given in Table 1. For each map, the first contour and the contour interval are the same. In the 1.2mm maps of L1498, 
L1495, L1400K, and L1517B, contours start at 5 mjy/ll"-beam, although the maps have been convolved to a resolution of 20”, and in the 
L1544 map, contours start at 20 mjy/13”-beam. The line maps represent integrated intensities including all hyperfine components. The first 
contour of each map ordered from left to right and from top to bottom are as follows (all in K km s“^, in the main beam temperature scale): 
L1498 (0.2, 0.075, 0.15, 0.3, 1.5); L1495 (0.3, 0.15, 0.45, 1.5); L1400K (0.2, 0.1, 0.3, 1.0); L1517B (0.2, 0.065, 0.1, 0.3, 1.5), L1544 (0.3, 0.13, 
0.15, 0.55, 2.0). Note that the point source to the NE of L1495 in the 1.2mm continuum map is IRAS 04112+2803 (aka CW Tau). 
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ers (dust plus the low dipole moment isotopomers C ^®0 
and C^^O), and the three bottom panels display maps 
of high density tracers (the large dipole moment species 
CS, N 2 H+, and NH 3 ). (The L1544 1.2mm map is from 
Ward-Thompson et al. 1999.) As it can be seen by 
simple inspection, there is a systematic pattern of emis¬ 
sion shared by all cores. On the one hand, the continuum, 
N 2 H+, and NH 3 maps are centrally concentrated, have ap¬ 
proximately the same peak position and share very similar 
shapes (nearly round in L1498, L1495, and L1517B, and 
more elongated in L1400K and L1544). The C^®0, C^^O, 
and CS maps, on the other hand, are much more diffuse 
and fragmented, and have maxima which do not coincide 
with those of the dust, N 2 H+, and NH 3 . In fact, many 
maps of the CO isotopomers and CS have relative minima 
at the position where the other molecules peak, suggest¬ 
ing almost an anti correlation between the two groups of 
tracers. 

The contrast between the centrally peaked dust emission 
and the fragmented C^®0 and C^^O maps (top rows in Fig. 
1 ) is especially striking, given that dust and C^® 0 /C ^’^0 
are both expected to trace the column density of the invis¬ 
ible H 2 component. We can rule out a major distortion in 
the CO isotopomer maps due to optical depth or satura¬ 
tion effects, since there is excellent agreement between the 
C ^®0 map and that of the rarer C^^O in each core where 
we have observed C^^O (Fig. 1). In addition, the intensity 
ratios between C ^®0 and C^^O are close to the expected 
isotopic ratio of 3.65, and a hyperfine analysis of the C^^O 
multiplet shows a negligible optical depth (< 0.1). We 
therefore conclude that both species are optically thin (see 
section 5.3 for a quantitative analysis). 

Temperature gradients cannot be causing the morpho¬ 
logical differences either, as the relative intensities of 
NHs)!,!) and (2,2) indicate a constant gas temperature of 
about 10 K in all observed sources (section 5.5). The cores 
moreover, being starless, lack internal heating sources to 
warm up the dust. Thus, we conclude that the difference 
between the continuum and C ^®0 maps has arisen from 
spatial differences in the dust and C ^®0 column densities. 
The excellent agreement between the dust emission and 
the emission from the high density gas tracers N 2 H+ and 
NH 3 (Fig. 1) suggests that the 1.2mm continuum is a 
faithful tracer of the gas component, especially at the core 
center, and that the CO isotopomers miss the high density 
central region. This must be the result of a drop in CO 
abundance at high density, similar to that found in sev¬ 
eral starless cores by other authors (L1498: Willacy et al. 
1998, IC5146: Kramer et al. 1999; Bergin et al. 2001, 
L977: Alves et al. 1999, L1544: Caselli et al. 1999, 
L1689B: Jessop & Ward-Thompson 2001) 

Before estimating the abundance profile of each molec¬ 
ular species using a full radiative transfer analysis, we 
present a simple argument to estimate a lower limit to 
the CO abundance drop at some of the core centers. For 
this, we use the maps in Fig. 1, which show that the 
C ^®0 and C ^’^0 emissions not only avoid the dust peak, 
but in most cases present a relative minimum at that po¬ 
sition. Simple radiative transfer shows that the intensity 
of an optically thin, thermalized line is proportional to 
the molecular column density, so at a point in the map 
with impact parameter 5, the CO isotopomer emission is 


proportional to 

J nHAr)^{r)dz, 

where nH 2 is the H 2 number density, X is the molecu¬ 
lar abundance, z is the length along the line of sight, and 
r = + z^. In order to reproduce the observed C^^O 

and C ^®0 emission minima toward the core peaks, the 
amount inside the integral has to reach a relative mini¬ 
mum, and this can only happen if the molecular abun¬ 
dance decreases to more than compensate the central den¬ 
sity peak. In other words, the observed C^^O and C^®0 
emission minima imply that the fraction of C ^®0 and C^^O 
decreases faster than l/n(r) at the core centers. 

4. CORE MODELING: METHOD DESCRIPTION AND 
DENSITY PROFILES 

4.1. Method Description and Justification 

In the following sections we present the radiative trans¬ 
fer modeling needed to derive the full chemical compo¬ 
sition of the cores. Before discussing all its details, we 
review here the main steps and assumptions involved in 
the process. 

A basic simplification in our calculations is that we 
model the cores as spherical systems. The cores were se¬ 
lected from previous work (Caselli et al. 2001a) for hav¬ 
ing nearly circular contour maps. The maps of continuum, 
N 2 H+, and NH 3 emission in Fig. 1 show that this is in fact 
the case in most systems (L1400K being the most deviant 
core). Although some cores would seem better modeled 
as spheroids, our lack of information on the line-of-sight 
component makes any guess about this third dimension 
highly uncertain. It seems, therefore, more convenient to 
use spherical models for the cores and to compare the pre¬ 
dictions of those models with azimuthal averages derived 
from the data. This will be our approach in what follows. 

The first step in our modeling is to derive density and 
temperature distributions for each core. We determine the 
density using the 1 . 2 mm continuum emission, which seems 
a reliable tracer (e.g., Andre et al. 1996; Ward-Thompson 
et al. 1999) and for which we have the data of highest an¬ 
gular resolution. To measure the gas temperature, we have 
carried out a standard LTE analysis of the NH 3 data and 
derived for each core a constant temperature of approxi¬ 
mately 9.5 K, except for LI544 for which we derive 8.75 
K. These values are in agreement with other estimates of 
these and similar cores (e.g.. Fuller & Myers 1993), and 
they will be later confirmed by our full, non LTE NH 3 
analysis in section 5.5. 

Once densities and temperatures have been estimated, 
we derive molecular abundance profiles by solving the line 
radiative transfer with a slightly modified version of the 
Monte Carlo code from Bernes (1979). For each core we 
vary the molecular abundance of each species and the gas 
velocity field until the model emission (convolved with the 
appropriate Gaussian beam) fits both the observed radial 
distribution of integrated intensity and the shape of the 
emerging spectrum toward the core center. This approach 
is as close as we can get from modeling the emission at each 
core position with our spherical models (see, e.g., Bieging 
& Tafalla 1993, for a similar procedure). 

As a result of the above modeling, we derive a set of 
physical and chemical parameters that, for each core, si- 
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multaneously fit all the data in Fig. 1 plus the emergent 
central spectra. Fitting such a large number of observa¬ 
tional data simultaneously imposes strict constrains in the 
model parameters, and this has made their determination 
a highly non linear process. Although the description in 
the following sections is sequential, the real fitting process 
required a number of iterations and back-and-forth mod¬ 
eling until self-consistency was achieved. 

4.2. Density Profiles 

We start our analysis by deriving a density distribution 
for each core, and we do this by modeling its dust 1 . 2 mm 
emission. As our models assume spherical symmetry, we 
first create spherical equivalents of each core by taking 
azimuthal averages of the data. For the roundest cores 
L1517B and L1495, we simply average the data in circles, 
and for the more elongated L1498, L1400K, and L1544, we 
average the emission along ellipses with aspect ratio (b/a) 
and position angle (PA) as indicated in Table 2. In these 
latter cases, the radial coordinate is the geometric mean 
of the major and minor axes of the averaging ellipse. The 
resulting profiles, both in linear and logarithmic scale, are 
presented in Figure 2. 

The goal of our modeling is to find a set of density distri¬ 
butions whose 1 . 2 mm dust emission fits the radial profiles 
shown in Fig. 2, and this requires first choosing a set of 
physical parameters for the dust grains. Following Andre 
et al. (1996), we use a 1.2mm emissivity Ki. 2 mm = 0.005 
cm^ g“^, although it should be noted that this parameter 
has an uncertainty of about a factor of 2. For dust temper¬ 
ature, we use a constant value equal to the gas temperature 
estimated from the NH 3 analysis (Sec. 5.5). Recent work 
by Evans et al. (2001) and Zucconi et al. (2001) suggests 
that the dust temperature in a core may decrease toward 
the center due to extinction in the warming interstellar 
radiation field. Our NH 3 data, however, suggests that the 
gas temperature is constant across the core (Sec. 5.5), and 
given that the central densities we derive (~ 10 ® cm“®) 
are high enough for gas and dust to be thermally coupled 
(Goldsmith & Langer 1978; Goldsmith 2000), we con¬ 
clude that if central cooling is present, it occurs at scales 
smaller than the resolution of our molecular line data, so 
it can be safely ignored in our models. 

If the dust kinetic temperature and the dust emissivity 
are constant, the emitted flux density from the core at the 
bolometer central frequency (240 GHz) is simply given by 

Sl,2mm — ^beam ^1.2mm 'bbl A^(i72) PuifPd) 

= 5.2 X 10“^^ B^{Td) [e{f')f N{H 2 ) mJy/beam, 

where flbeam is the telescope beam solid angle, m is the 
mean molecular mass, N{H 2 ) is the H 2 column density 
(in cm“^ in the lower equation), 6{") is the FWHM of 
the beam in arcsec, and B^{Td) is the Planck function at 
the dust temperature T^. With this equation, finding the 
core density law becomes simplified to finding the func¬ 
tion n{r) that gives rise to the appropriate N{H 2 ) profile. 
Unfortunately, the process is complicated by the need to 
take into account the beam smoothing and spatial filtering 
introduced by the on-the-fly bolometer observation. This 
problem, discussed in detail by Motte & Andre (2001), 
requires that a simulation of the on-the-fly observation is 
carried out before comparison with the data. In our case, 


we have used the NIG software, which lets us simulate an 
on-the-fly observation of the model with the same param¬ 
eters as those used during the telescope run, and reduce 
the simulated data the same way the real data were re¬ 
duced. This guarantees that data and model are properly 
compared. 

To find the density distributions that reproduce the ob¬ 
served radial profiles, we iterated the fitting procedures 
until reasonable convergence. Previous work has shown 
that single power-law density profiles do not fit the emis¬ 
sion from starless cores, and that a central flattening is 
always needed to reproduce the data (Ward-Thompson et 
al. 1994; Andre et al. 1996; Bacmann et al. 2000; Alves 
et al. 2001). Thus, it has become a standard practice to 
use double power laws with the inner portion nearly flat. 
The discontinuous derivative of these fits, however, seems 
rather artificial, and it appears more likely that real cores 
have smoother density profiles, with the double power-law 
expressions being just an approximation. For this reason, 
we have searched for a family of analytic density profiles 
that combine the power-law behavior for large r and the 
central flattening at small r. After different tests, we have 
chosen profiles of the form 

/ \ ^0 

"(’■) = 1+ 

where no is the central density, ro is the radius of the “flat” 
region ( 2 ro is the full width at half maximum), and a is 
the asymptotic power index. As Fig. 2 shows, this family 
of curves, with the values for no, ro, and a indicated in 
Table 2, fits each of the observed radial brightness distri¬ 
butions with great accuracy, and it therefore constitutes 
the basis of our core modeling. Values of no range be¬ 
tween about 10 ® and 10 ® cm“®, and their major source of 
uncertainty is their linear dependence on Ki, 2 mm (which 
has a factor of 2 uncertainty), tq is probably the best de¬ 
termined parameter and ranges from 3,000 to 10,000 AU, 
with an uncertainty of about 1,000 AU. Finally, the a pa¬ 
rameter ranges from 2 to 4, and is probably accurate to 
0.25-0.5. (Whitworth & Ward-Thompson (2001) have re¬ 
cently proposed a similar expression, although they require 
a > 4, which seems outside the range of our empirically 
determined values. See also Langer & Willacy (2001) for 
an alternative fit to the L1498 density profile.) 

5. ABUNDANCE PROFILES: EVIDENCE FOR CHEMICAL 
DIFFERENTIATION 

5.1. Monte Carlo Model Parameters 

Once the core density profiles have been determined 
from the continuum emission, we derive molecular abun¬ 
dances for all the observed species by fitting their distri¬ 
bution of line intensities. We do this by solving the equa¬ 
tions of radiative transfer with a Monte Garlo code (Bernes 
1979), and as mentioned before, we assume a gas kinetic 
temperature of 8.75 K for L1544 and 9.5 K for the rest of 
the cores (section 5.5 shows how these temperatures are 
required by the NH 3 data). 

For each core, the free parameters in our Monte Garlo 
models are the molecular abundances, the gas velocity field 
(both systematic and turbulent), and the core maximum 
radius (almost immaterial given the steep density profiles). 
The velocity field is well constrained by the need to fit si- 
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Fig. 2.— Radial profiles of 1.2mm continuum emission in both linear and log-log scale. (L1544 data from Ward-Thompson et al. 1999.) 
Solid squares are telescope measurements and lines are models for which a simulation of an on-the-fly observation has been done (see text). 
For L1495 and L1517B, the profiles represent circular averages, while for L1498, L1400K, and L1544, they represent elliptical averages with the 
parameters shown in Table 2. See also Table 2 for the input parameters of the models. Scale is in mjy/ll"-beam for L1498, L1495, L1400K, 
and L1517B, and mjy/13^^-beam for L1544. In contrast to Fig. 1, the data from the first four cores have not been convolved additionally 
after the telescope observation. 
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multaneously the central line shape of C^®O(l-0), CS(2- 
1), N 2 H+( 1 - 0 ), NH 3 ( 1 , 1 ), and NH 3 ( 2 , 2 ). In this way we 
find that all cores but L1544 are consistent with a con¬ 
stant turbulent linewidth, with a value of FWHM = 0.167 
km s“^ (variable vturb = 0.1 km s“^ in the Bernes 1979 
Monte Carlo program). For L1544, we model the emission 
with a turbulent linewidth of 0.13 km s“^ at the outermost 
radius increasing inward as up to a value of 0.25 km 
s“^ (Tafalla et al. 1998 also found an inward increase of 
the turbulent linewidth, and recently Caselli et al. 2001c 
have found a similar pattern studying the linewdith in the 
plane of the sky). This larger value of the linewidth in 
L1544 probably reflects the complicated kinematics of this 
core, already discussed by Tafalla et al. (1998) and Caselli 
et al. ( 2001 c). 

Together with the turbulence velocity field, we And that 
some cores have systematic velocity gradients along the 
line of sight. These gradients are required to fit the dis¬ 
placement of the CS self absorptions with respect to the 
systemic velocity (section 5.2), and, for simplicity, we 
parametrize the velocity held with a linear function. The 
values we find for dV/dr (Table 3) are consistent with the 
typical core velocity gradients across the line of sight mea¬ 
sured by other authors (e.g., Goodman et al. 1993). 

The maximum core radius, as mentioned before, is not 
a critical parameter in our modeling, because the steep 
density gradient of each core sets an almost natural outer 
edge. All cores but L1544 are well fit with a maximum 
core radius of 190' (4 x 10^^ cm at the distance of Taurus 
and 5 x 10^"^ cm for L1400K), and L1544 requires twice 
that value. This large radius for L1544 is needed in or¬ 
der to fit its very deep CS self absorption, a result already 
found by Tafalla et al. (1998). These authors carried out 
a similar Monte Carlo fitting of the CS lines in L1544, but 
ignored possible abundance variations with radius. The 
models presented in this paper supersede those in Tafalla 
et al. (1998). 

A number of tests have been carried out to ensure that 
the Monte Carlo results presented here are independent of 
the internal program parameters. During searches for the 
best fits, the core was divided into 100 shells spaced loga¬ 
rithmically, the number of model photons was set to 2000 , 
40 iterations were performed, and no reference radiation 
field was used. Previous experience had shown that these 
parameters were sufficient to avoid dependence on the ini¬ 
tial conditions and to reach a stable solution. As a further 
test, all best models for L1517B (one of the most regular 
cores) have been repeated doubling the number of shells 
and iterations and increasing the number of photons by 
an order of magnitude; no significant difference has been 
found between these new runs and the standard cases. We 
take this stability of the results as an indication that the 
outputs from our models are true solutions to the molecu¬ 
lar radiative transfer problem, and are not affected by the 
internal details of the Monte Carlo code. 

5.2. and CS Abundances 

C^®0 and CS are closed-shell linear rotors with known 
collisional rates, and this makes their radiative transfer 
the easiest to calculate among all our species. For C^®0, 
we use the molecular constants from Lovas & Krupenie 
(1974) together with the rate coefficients of collisional ex¬ 


citation by para H 2 from Flower & Launay (1985). For 
CS, we also use the molecular constants from Lovas & 
Krupenie (1974), together with the collisional coefficients 
for para H 2 from Green & Chapman (1978). A total of 
9 rotational levels are considered for each molecule, and 
the resulting model line intensities are convolved with a 
Gaussian beam to simulate the result of an observation 
with the FCRAO telescope (FWHM of 46" for C^® 0(1-0) 
and 50" for CS(2-1), according to Ladd & Heyer 1996). 

Our first set of models assumes spatially constant abun¬ 
dances for both C^®0 and CS. For C^®0, we assume an 
abundance of 1.3 x 10“^ for all cores and for CS we use dif¬ 
ferent abundances so that we approximately fit the emis¬ 
sion at large radius (> 100") for each core. The resulting 
radial profiles, indicated by dashed lines in the left panels 
of Fig. 3, fail to reproduce the observed intensities near 
the core center by more than a factor of 2 and are, there¬ 
fore, inconsistent with the data. A similar set of constant 
abundance models with an abundance low enough to fit 
the central emission (models not shown in Fig. 3) also 
clearly fails to fit the C^®0 and CS radial profiles, in this 
case by underestimating the emission at large radius. 

In order to improve the fits, we have run additional mod¬ 
els with a central abundance drop. In section 3 we have 
seen that the C ^®0 abundance has to decrease faster than 
l/n(r), so with the Monte Carlo program, we have ex¬ 
plored different analytical expressions of increasing steep¬ 
ness: l/n(r)^, an exponential form, and a central abun¬ 
dance hole. The results derived from these models are 
practicably indistinguishable due to the limited resolution 
of the FCRAO data (~ 50"), so, for simplicity, we will 
only present here the results derived with the exponential 
abundance law: 

X{r) = Aoexp(-n(r)/nd), 

where the free parameters Xq and Ud represent the low- 
density abundance limit and the e-folding size of the abun¬ 
dance drop. A more detailed discussion of the shape of the 
abundance profiles using higher angular resolution data 
and multiple transitions will be presented in a forthcom¬ 
ing paper (Tafalla et al. 2001). 

The results from this last set of models are indicated 
by solid lines in Fig. 3. For C^®0, all models have the 
same Xq value of 1.7 x 10“^, following Frerking et al. 
(1982), and the Ud parameter ranges from about 1.5 x 10^ 
to 5.5 X lO"* cm“®, depending on the core (see Table 4). 
For CS, different cores need slightly different values of Xq, 
but most estimates are close to 4 x 10“® (Table 4). The 
Ud values for CS are also of the order of a few 10"^ cm“® 
(10® cm“® in L1544). 

As seen in Fig. 3, the exponential-law abundance mod¬ 
els fit the observed radial profiles better than the constant 
abundance models, especially in the inner 100". For L1498 
and L1517B, the fits are good at all radii, and the same 
could be said for L1544, but the scatter in the radial pro¬ 
files makes it difficult to choose the best fit at large radius. 
For L1495 and L1400K, however, the fits deviate from the 
data, and the scatter of the radial profiles makes it impos¬ 
sible to find a perfect fit even if we could draw an arbitrary 
radial profile. The origin of these fitting problems can be 
understood by looking at the maps in Fig. 1, which show 
that, in these two cores, the C®®0 and CS emission is 
so different from that of the dust, NH 3 , and N 2 H+ that it 
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Fig. 3a. — Radial profiles of integrated emission (left) and central emerging spectrum (right) for C^®O(l-0) and CS(2-1) in each of our 
cores. Observations are indicated by solid squares in the radial profiles and by histograms in the spectra (all temperature units are in the 
main beam scale). The results from two radiative transfer models are plotted over the radial profiles: the dashed lines are constant abundance 
models and the solid lines are models with a rapid decrease in the central abundance. (Constant abundance models assume X(C^^O) = 1.3 
10“^ in all cores and X(CS) = 3.0 10“® in L1495 and L1517B, 2.0 10“® in L1498 and L1544, and 6.0 10“® in L1400K. See Table 4 for 
parameters of the models with a central abundance drop.) Note how the constant abundance models fail to fit the central emission, while 
the models with variable abundance fit the data at all radii satisfactorily. The predicted emerging intensities of these variable abundance 
models (after convolution) are presented by solid lines overlaid on the central spectra in the right panels. Note the presence of additional 
velocity components in the emission toward L1498 and in and CS toward L1495. These components probably arise from unrelated 

ambient gas and have not been included in the integrated emission of the radial profiles. 
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seems to trace additional structures unrelated to the cores. 
These structures most likely arise from nearby lower den¬ 
sity gas, and their modeling is outside the reach of our 
simple scheme, which uses the continuum emission as the 
basis of the core model and assumes spherical symmetry. 
This last assumption, in addition, has the most extreme 
deviation in the case of L1400K, which is clearly elongated 
NE-SW and has additional extensition to the west (Fig. 
1). The above limitations in the L1495 and L1400K mod¬ 
els, however, do not contradict the idea of a C^®0 and CS 
abundance decrease in these two cores. Indeed, it seems 
that their central abundance decrease is more severe than 
that shown by other cores, and it is that very fact which 
makes them almost invisible in C^®0 and CS. 

An additional requirement for our models is that they 
also fit the central spectrum, both in integrated intensity 
and detailed shape. As the right panels of Fig. 3 show, 
this is the case for L1498, L1517B, and L1544, while there 
are differences between model and data in L1400K and 
L1495 (note that the extra wing in the L1498 C^®0 spec¬ 
trum most likely arises from the extra red component seen 
in the maps of Lemme et al. 1995). As with the radial 
profiles, the worst fits occured in the cases of L1495 and 
L1400K. This probably results from the difficulty in mod¬ 
eling extended, non spherically-symmetric material, a sit¬ 
uation further complicated in L1495 by the presence of gas 
at different velocities (see our C^®0 spectrum and Good¬ 
man et al. 1998 for a detailed discussion of the velocity 
structure in L1495). For the cores with good fits (L1498, 
L1517B, and L1544), the agreement between model and 
data is achieved by finding the appropriate combination of 
systematic and turbulent velocity fields, a process carried 
out simultaneously with the modeling of the N 2 H“'" and 
NH 3 emission (sections 5.4 and 5.5) because our model re¬ 
quires that, for each core, all lines be reproduced with a 
single velocity held. 

The CS spectra are particularly sensitive to the system¬ 
atic velocity field due to the presence of self absorptions. 
These features naturally arise from the combination of a 
drop in excitation caused by the density fall off with radius 
and a relatively large abundance of CS in the outer core 
layers. Regarding the modeling of the emission, the dif¬ 
ference in velocity between the self absorption dip and the 
line center constrains any possible velocity gradient along 
the line of sight. Most of our cores show evidence of such 
gradients, which we have modeled for simplicity using lin¬ 
ear functions with the coefficients given in Table 3. We find 
both inward and outward velocity fields among our sample 
cores, with the best examples of each case being L1498 (in¬ 
ward) and L1517b (outward). We note that our model for 
the infall case L1544 lacks a gradient along the line of sight. 
This occurs because our modeling process was restricted to 
fit only the shape of the central spectrum, which appears 
rather symmetric and deeply self-absorbed (see Fig. 3). 
As Tafalla et al. (1998) have shown, L1544 presents ev¬ 
idence of extended inward motions in CS when off-center 
positions are considered, an effect that we cannot model 
here because of our assumption of spherical symmetry. 

To conclude our analysis of the C^®0 and CS emission, 
we will estimate the abundance contrast between the core 
center and the outer layers implied by our models. To do 
this, we define for each core the parameter / as the ra¬ 


tio between the model abundances at r = 20 " (half our 
beam) and r = 100". This / parameter, shown in Table 
4, is typically of the order of 10“^ for both C^®0 and CS, 
indicating abundance drops for these two molecules of two 
orders of magnitude in each core. Such strong abundance 
drops are, in fact, consistent with a case in which the cores 
have no C^®0 and CS molecules at their centers. 

5.3. C^'^0 and Data 

The Monte Carlo modeling takes into account optical 
depth and trapping effects, so the poor fit of the constant 
abundance models in Fig. 3 cannot be explained by a fail¬ 
ure to take those factors into accout. Still, the C^®O(l-0) 
lines have moderate optical depth and the CS(2-1) spectra 
are obviously thick and self absorbed, so it seems appro¬ 
priate to test the results from the previous section using 
optically thin lines of C^^O and C^^S. Here we carry out 
such testing with C^^O(l-O) observations of the cores with 
brightest C^®0 emission (L1498, L1517B, and L1544, see 
Fig. 1) and C^^S(2-1) data from the most self absorbed 
core (L1544). 

The C^^O molecule has hyperfine structure, and this 
in principle could complicate its radiative transfer mod¬ 
eling. The low optical depth and dipole moment of this 
molecule, however, causes its level excitation to be dom¬ 
inated by collisions, leaving trapping and splitting effects 
negligible. This allows us to model the C^^O radiative 
transfer with the Monte Carlo code as if there were no 
hyperfine structure, and to predict the combined inten¬ 
sity of the three hyperfine components. The lack of model 
prediction for the individual components is of little conse¬ 
quence, as the low signal-to-noise ratio of our C^^O data 
makes it necessary to combine the observed components 
before comparing them with the model predictions (see 
sections 5.4 and 5.5 for model line predictions for species 
with hyperfine structure like N 2 H+ and NH 3 ). 

To test the C^®0 calculations from the previous section, 
we need to create equivalent models for the rarer C^^O. As 
the relative abundance of C ^®0 over C^^O seems constant 
over the Galaxy with a well measured value of 3.65 (Pen- 
zias 1981), the C^®0 abundance curves determined pre¬ 
viously can be automatically converted into C^^O abun¬ 
dance profiles just by dividing by 3.65. This means that 
our C^^O Monte Carlo modeling is a zero-free-parameter 
calculation, as it has no internal values that can be used 
to improve the fitting. 

The top panels of Fig. 4 present a comparison between 
the results of our C^^O(l-O) models and the FCRAO ob¬ 
servations. Here, once again, the dashed lines represent 
constant abundance models and the solid lines indicate 
models with an exponential abundance drop. As before, 
the models with constant abundance fit only the outer core 
emission and deviate dramatically from the observed in¬ 
tensity toward the center. The discrepancy between the 
C^^O constant abundance models and data is larger than 
seen for C^® 0 , due to the lower optical depth of the rarer 
isotope, and reaches values of up to a factor of 5 toward 
the core centers. 

In contrast to the poor fit of the constant abundance 
models, a remarkable fit is obtained with the models that 
have an exponential abundance drop (Fig. 4). These mod¬ 
els have the same abundance profiles as those used for 
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Fig. 4.— Radial profiles of integrated emission for C^^O(l—0) (three top panels) and C^‘^S(2-1) (bottom panel) toward selected cores. 
Observations are indicated by solid squares and results from radiative transfer models are indicated by lines. The models assume the same 
abundance profiles as those in Fig. 3, with the ISM isotopic ratio of 3.65 for C^^O/ (Penzias 1981) and the Solar ratio of 22.7 for 

C^^S/C^'^S. Note how the constant abundance models (dashed lines) miss the data toward the core centers by large factors. The models with 
central abundance drops (solid lines), in contrast, automatically fit the data without adjusting any additional parameter. 
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CISq (Table 4), with the only difference being that the 
Xq value has been divided by 3.65. The good match of 
these models, both at large and small radii, constitutes 
a further and final proof that the abundance of the CO 
isotopes decreases sharply and systematically at high den¬ 
sities in all cores of our sample. 

Our modeling of the C^‘*S(2-1) emission is also a zero- 
free-parameter calculation. The C^^S/C^'^S isotopic ratio 
in the interstellar medium equals, to a good approxima¬ 
tion, the Solar value of 22.7 (Lucas & Liszt 1998), so the 
C^'^S abundance profiles are just 22.7 times smaller than 
the CS abundance profiles. To test the main isotope CS 
calculations in LI544, we now run Monte Carlo radiative 
transfer models for C^^S(2-1) for the cases of constant 
abundance and exponential drop. 

As Fig. 4 shows, the model with constant abundance 
(dashed lines) is again inconsistent with the data, as it 
overpredicts the central intensity by more than a factor of 
3. The model with an exponential abundance drop, on the 
other hand, again fits the data remarkably well, predict¬ 
ing a radial profile of the correct intensity and shape. This 
confirms our conclusion that the CS abundance decreases 
rapidly toward the core center and shows that our Monte 
Carlo modeling can predict accurate abundance profiles 
even in a case with large optical depth and strong self 
absorption. 

5.4. N 2 H+ Data 

Solving the radiative transfer for N 2 H''' is much more 
complicated than for CO and CS, as the two N atoms of 
N 2 H+ make its rotational level structure split into multi¬ 
ple hyperfine components (7 for J=l-0, 38 for 2-1, 45 for 
3-2, etc.). The large number of resulting transitions, and 
their possible overlap, makes the solution of the radiative 
transfer a formidable task, and places it outside the scope 
of our present analysis. In this section, we discuss how to 
simplify the N 2 H“'' radiative transfer treatment, and how 
to derive reasonably-accurate molecular abundance pro¬ 
files for all our cores. 

The presence of hyperfine (hf) splitting introduces two 
new elements into the radiative transfer. First, the split¬ 
ting allows the relative population of the hyperfine sub- 
levels to depart from their statistical weights. This effect, 
which gives rise to hyperfine multiplets with “non-LTE” 
intensity ratios, has been observed in N 2 H+(l- 0 ) toward 
starless cores by Caselli et al. (1995). Although well- 
detected, this effect involves less than 10% of the emit¬ 
ted 1-0 radiation (Caselli et al. 1995), so it must repre¬ 
sent only a small perturbation in the total emerging flux. 
In the simplified treatment presented here, we will ignore 
this effect assuming that inside each rotational transition 
the sublevels are populated according to their statistical 
weights. 

The second effect of the hyperfine splitting is a change 
in the line trapping. Because of the splitting, the pho¬ 
tons from each rotational transition can escape more eas¬ 
ily, and this decreases their role in excitation. This effect 
will only be significant if trapping dominates collisions as 
the main excitation mechanism, and this can only occur 
at very large optical depths. A simple LTE fit to the ob¬ 
served (1-0) multiple! (the one expected to be thickest) 
shows that, in our cores, the averaged optical depth of the 


hf components does not exceed 2 even at the core maxi¬ 
mum, and that in most cases the lines are optically thin. 
This suggests that collisions dominate the excitation, and 
that the trapping decrease represents only a minor contri¬ 
bution. To further explore this point, we have run radia¬ 
tive transfer models using two extreme limits. In one, we 
completely ignore the hf splitting, so the trapping is maxi¬ 
mum, and in the other, we artificially turn off the trapping 
so it has no effect on excitation. The real case will nec¬ 
essarily lie between these two limits, closer to one or the 
other depending on optical depth. As a result of this test, 
we find that for the conditions of our cores, both limits 
predict the same shape of the radial profile of integrated 
intensity. Thus, the shape of the N 2 H“'' abundance curve 
we derive from fitting the radial profile is independent of 
our treatment of the hf splitting. The absolute value of 
the abundance, in addition, typically differs between the 
two cases by just 30% (80% in L1495), indicating that 
the effect of the treatment of the hf splitting is small, and 
probably comparable to the uncertainties in the collisional 
coefficients (see below). 

Given the small effect of the hf splitting in the exci¬ 
tation, we can simplify the radiative transfer calculation 
by separating the problem in two steps: the solution of 
the level excitation and the prediction of the emergent 
spectrum. In the first step, we use the Monte Carlo code 
without hf splitting, and calculate the combined popula¬ 
tion of each rotational level. Eor the collisional rate coeffi¬ 
cients, we use the recent set for collisions between HCO’*' 
and para H 2 from Elower (1999). As Monteiro (1984) 
has shown, the similar rotational constants and electronic 
structures of N 2 H+ and HCO“'" imply that the two species 
have similar cross sections, and they do in fact have colli¬ 
sional rates with He which are indistinguishable (see also 
Green 1975). Until a specific calculation for collisions 
between N 2 H+ and H 2 exist, the use of the HCO"'' rates 
remains the most accurate choice. 

Once the combined population of the rotational levels 
has been derived with the Monte Garlo code, we take 
into account the full hf structure to integrate the equa¬ 
tion of radiative transfer along the line of sight and predict 
the emergent 1-0 spectrum. We assume that the hf sub- 
levels are populated according to their statistical weights 
and use the relative line frequencies from Gaselli et al. 
(1995). This emerging spectrum is finally convolved with 
a Gaussian function of 52.5 arcsec FWHM, to simulate the 
smoothing effect of the FGRAO telescope beam. 

The results of the N 2 H+ calculations are presented in 
Fig. 5, where constant abundance models are compared 
with FGRAO observations for each core in our sample. As 
Fig. 5 shows, models with constant N 2 H“'' abundance pro¬ 
vide satisfactory fits in most cases to both the radial profile 
of integrated intensity (left panels) and the central emerg¬ 
ing spectrum (right panels), in sharp contrast to what was 
found for G^®0 and GS. Although there are small devia¬ 
tions between model and data in Fig. 5, they most likely 
result from our simplified modeling and do not represent 
real abundance variations across the cores. This seems to 
be the case in L1400K, where the bright emission at large 
radius (and overal larger scatter) is due to the elongated 
shape of this core and to the presence of an extra N 2 H“'' 
component to the west (see Fig. 1). Also, our inability to 
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Fig. 5.— Radial profiles of integrated emission (left) and central emerging spectrum (right) for N2H“^(1—0). Observations are indicated 
by solid squares in the radial profiles (sum of all hyperfine components) and by histograms in the spectra (all temperature units are in the 
main beam scale). Constant abundance models are indicated by solid lines, and unlike the case of and CS (Fig. 3), they accurately fit 

the data at all radii (see parameters in Table 4). For the emerging spectra, only the central hyperfine group is presented to show in detail 
the good matching of both the linewidth and the hyperfine ratios. To use the same velocity scale in all the spectra, both data and model for 
L1400K have been shifted by 4 km s~^, and by 1 km s~^ for L1517B. 
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fit the non-gaussian shape of the L1544 spectrum is due 
to the presence of two velocity components in this core 
(Tafalla et al. 1998; Caselli et al. 2001c), which are 
simply modeled here with a broader gaussian component. 
Taking into account these effects (and the scatter in the 
radial profiles), we conclude that our data are consistent 
with all cores having constant N 2 H''' abundance, which we 
estimate to be of the order of 10“^° (see Table 5 for the 
value derived for each core). 

Although proving that the N 2 H+ abundance is constant 
has required a full radiative transfer calculation, it is clear 
from Fig. 1 that the similarity between the 1.2mm contin¬ 
uum and the N 2 H~'' data implies an abundance behavior 
of this molecule different from that of CO and CS. Caselli 
et al. (1999) and Bergin et al. (2001) have reached simi¬ 
lar conclusions for the cases of L1544 and 1C5146, respec¬ 
tively, using a simpler (standard LTE) analysis than the 
one presented here. The different behavior of N 2 H+, on 
the one hand, and C^®0 and CS, on the other, is the first 
element of a chemical differentiation pattern that affects 
all our cores, and a clear sign that the common assump¬ 
tion that different molecular species trace the core interior 
equally needs serious revision. 

5.5. NH 3 Data 

The NH 3 molecule is a symmetric top with inversion 
doubling and consists of ortho and para species, which co¬ 
exist almost independently of each other (see, e.g., Ho & 
Townes 1983). The (J,K) = (1,1) and (2,2) inversion lines 
we have observed arise from para-NHa, so our modeling 
will only deal with this species. 

Like N 2 H+, NH 3 has hyperffne splitting, thus, our treat¬ 
ment of this effect again involves some simplifications. For¬ 
tunately, the effect of the hf splitting in NH 3 is even less 
critical than in N 2 H+. On the one hand, NH 3 does not 
present hf anomalies toward dark clouds (Stutzki et al. 
1984, also our own data), so the population of the hf sub- 
levels should be proportional to their statistical weights. 
On the other hand, the decrease of trapping due to the 
hf splitting has very little effect on the collisionally dom¬ 
inated NH 3 excitation. This once again, can be seen by 
comparing models with the two extreme trapping treat¬ 
ments (see section 5.4), which for L1517B (the roundest 
and second brightest core in NH 3 ) give rise to an abun¬ 
dance difference of less than 10%. As before, such a differ¬ 
ence is a negligible factor, comparable to our calibration 
uncertainty. 

Following our N 2 H+ treatment, we first solve the NH 3 
level populations assuming no hf splitting and then cal¬ 
culate the emerging spectrum using sublevel populations 
proportional to the statistical weights. To do this, we 
have modified the Monte Carlo code to include the inver¬ 
sion doubling of the rotational levels in a symmetric top 
molecule. Given the low gas temperature in our cores (9.0- 
9.5 K), more than 99.99% of the para-NH 3 population is 
expected to lie in the 6 levels resulting from the inversion 
of the rotational states (J,K)=(1,1), (2,1), and (2,2), so 
only those will be considered in our models. The energies 
of these levels have been calculated with the parameters 
of Poynter & Kakar (1975), and the results have been 
checked against the values in Pickett et al. (1998). The 
6 inversion levels are connected by 5 radiative transitions. 


and their Einstein coefficients are given by the standard 
expression for a symmetric top molecule (see, e.g., Townes 
& Schawlow 1955) with a dipole moment of 1.476 debye 
(Cohen & Poynter 1974). Collisional coefficients connect¬ 
ing all levels have been calculated by Danby et al. (1988) 
for temperatures between 15 and 300 K, and their results 
show that the de-excitation coefficients depend weakly on 
temperature. Thus, we have used for our modeling the 15 
K de-excitation coefficients with para H 2 , and calculated 
the temperature-dependent excitation rates using the con¬ 
dition of detailed balance for the correct temperature. As a 
final step in our calculation, the emergent intensities have 
been convolved with a Gaussian of 40" FWHM, to simu¬ 
late the result of an observation with the 100 m telescope. 

As mentioned before, a simple LTE analysis of the NH 3 
data indicates a constant temperature for all cores. With 
our Monte Carlo calculation we can now test this result 
by fitting the ( 1 , 1 ) line and checking whether the ( 2 , 2 ) 
line fits automatically. Given that the ratio between the 
intensity of these lines is very sensitive to the kinetic tem¬ 
perature (Walmsley & Ungerechts 1983), any mismatch in 
the ( 2 , 2 ) model will indicate a deviation from the chosen 
temperature. 

The fit of the C^®O(l-0), CS(2-1), and N 2 H+(l- 0 ) spec¬ 
tra in the previous sections has already constrained the 
turbulent and systematic velocity fields in the models, and 
this leaves the molecular abundance as the only free pa¬ 
rameter in our NH 3 calculation. As we did for the other 
species, we begin the analysis by running models with con¬ 
stant NH 3 abundance and comparing them with the ob¬ 
served radial profiles. The dashed curves in Fig. 6 show 
that the emission from these models is too flat to repro¬ 
duce the observed central increase in both the ( 1 , 1 ) and 
(2,2) transitions, suggesting that constant NH 3 abundance 
models are inconsistent with the data (the L1400K emis¬ 
sion can be fit with a constant abundance model, but our 
data are too noisy to tell). 

The poor fit of the constant abundance models affects 
both the ( 1 , 1 ) and ( 2 , 2 ) transitions, so it cannot have re¬ 
sulted from a wrong choice in the gas kinetic temperature. 
It cannot be due either to an incorrect estimate of the ex¬ 
citation temperature, as the ( 1 , 1 ) and ( 2 , 2 ) lines are ther- 
malized to a good approximation (gas densities are close 
to 10 ® cm“® while the critical density is 2 x 10 ® cm“®). 
Different tests assuming constant excitation temperature 
or an excitation temperature as given by a two level ap¬ 
proximation (e.g., Stutzki & Winnewisser 1985) confirm 
the Monte Carlo results and show that fitting the observed 
radial profiles requires an enhancement of the central NH 3 
abundance. 

To model the NH 3 abundance enhancement we have 
tried different analytic expressions, such as power laws and 
exponential forms. The limited resolution of the observa¬ 
tions (40") makes it possible to fit the data with different 
expressions, as long as the abundance increases by a fac¬ 
tor of a few toward the center. Here we present the results 
of power law models (A(r) = Xq (n(r)/no)^), shown by 
solid lines in Fig. 6 . The values for Xq and l3 are given in 
Table 5, which shows that the abundance enhancement / 
ranges between 3 and 12 (excluding L1400K). 

The abundance enhancement models not only fit the 
( 1 , 1 ) radial profiles, but those of ( 2 , 2 ) as well, and this 
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Fig. 6 .— Radial profiles of NH3(1,1) a^nd (2,2) integrated emission (left) and central emerging NH 3 ( 1 , 1 ) and (2,2) spectra (right). Observa¬ 
tions are indicated by solid squares in the radial profiles (sum of all hyperfine components) and by histograms in the spectra (all temperature 
units are in the main beam scale). Dashed lines represent constant abundance models and solid lines represent models with a central abun¬ 
dance enhancement (see Table 5 for values). Note how in all cores but L1400K the models with constant NH 3 abundance fail to fit the central 
intensity (models have been set to fit outer core emission), while the models with a central abundance enhancement fit the data at all radii 
for both the (1,1) (2,2) transitions. For the emerging spectra, only part of the hyperfine group is presented in order to show in detail the 

matching of the linewidth and the hyperfine contrast. To use the same velocity scale in all spectra, both data and model for L1400K have 
been shifted by 4 km s“^, and by 1 km s“^ for L1517B. 
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implies a correct choice for the kinetic temperature of the 
gas from the LTE analysis. This is not surprising given 
the high degree of thermalization of the NH 3 levels, which 
makes the LTE approximation rather accurate. Although 
it is still possible that our 40"-resolution observations miss 
gas in the inner cores at significantly lower temperatures 
(Evans et al. 2001; Zucconi et al. 2001), the sharp cen¬ 
tral abundance increase, which biases the emission toward 
the innermost gas, suggests that this cold region has to 
be reduced to a radius significantly smaller than 0.015 pc 
(22" at 140 pc). 

Another satisfactory aspect of the NH 3 model is the 
good fit to the central spectrum, in particular, to its hf 
structure. This can be partially seen in the right panels of 
Fig. 6 , where the main and the first satellite components 
of the observed ( 1 , 1 ) spectra are compared with the model 
predictions. Given the absence of any free parameter to 
control the hf ratio of the model results, it is remarkable 
that in all cases the model predicts the correct intensity 
contrast between hf components, and therefore, the correct 
optical depth. This threefold fit of the NH 3 data (radial 
profiles of two transitions and hf ratio) strongly suggests 
that the abundance increase of NH 3 toward the core cen¬ 
ters is an inescapable consequence of the observations. 

6 . DISCUSSION 

6.1. Molecular Depletion and Comparison with Models 

The analysis of the previous sections shows that the 
cores in our sample share a similar pattern of chemical 
differentiation. All systems present a sharp, order of mag¬ 
nitude, central abundance drop in C^®0 and CS, a con¬ 
stant abundance of N 2 H+, and a central enhancement of 
NH 3 . This pattern seems independent of core properties 
like the gas velocity structure (inward/outward), the core 
size, or its exact shape, indicating that it has to arise from 
a robust and general chemical process. 

The most natural explanation for the abundance drops 
of C^®0 and CS is the depletion of these molecules onto 
cold dust grains at high densities. Molecular depletion 
has been widely expected and sought (e.g., Mundy & Mc- 
Mullin 1997), and our data add to an increasing number 
of observations that indicate its final detection (Kuiper et 
al. 1996; Willacy et al. 1998; Kramer et al. 1999; Alves 
et al. 1999; Caselli et al. 1999; Jessop & Ward-Thompson 
2001; Bergin et al. 2001). The systematic pattern in our 
sample cores shows that molecular depletion is not just 
limited to a few special cases, but that it characterizes 
the majority (or totality) of the starless core population 
at densities larger than a few 10^ cm“^ (Table 4). 

The coexistence of strong C^®0 and CS depletion with a 
constant N 2 H+ abundance may appear at first surprising, 
and even more so, the enhancement of NH 3 abundance at 
the density peak. A likely explanation for this behavior 
comes from the recent work by Bergin & Langer (1997, 
BL97 hereafter), who have presented a chemical model in 
which differential depletion arises naturally from the dif¬ 
ferent binding energies of the different species (see also 
Charnley 1997; Aikawa et al. 2001; Caselli et al. 2001c). 
According to BL97, the N 2 molecule binds more weakly 
to grains than do CO and CS, and thus it is more easily 
desorbed. This difference in binding energies gives rise to 
a pattern of depletion in which N-rich molecular species 


survive in the gas phase at higher densities than do CO 
and CS, which disappear at densities of about 10^ cm“^. 
To further explore the BL97 model predictions, we have 
used our radiative transfer calculation to estimate emerg¬ 
ing intensities for the cores in our sample. 

BL97 have studied the chemical evolution of a gas par¬ 
cel as its density increases with time following two dif¬ 
ferent models, one based on theoretical work by Basu & 
Mouschovias (1994, BM94 hereafter), and the other based 
on a phenomenological (accretion) fit to the contraction 
history of L1498 (Kuiper et al. 1996). The main differ¬ 
ence between these two models is that the BM94 contrac¬ 
tion is rather slow (as it is driven by ambipolar diffusion), 
while in the accretion fit model the density initially in¬ 
creases much faster. For both models, BL97 studied the 
effect of two types of grain mantles (CO and H 2 O), but a 
look at their abundance predictions shows that significant 
CO depletion only occurs in models having dust grains 
with (tightly bound) H 2 O mantles. Thus, in the following 
discussion we concentrate on models with these types of 
grains. 

To compare the BL97 predictions with our observations 
of real cores, we create a model core with the density dis¬ 
tribution of one of our observed cores and an abundance 
distribution that matches point by point the abundances 
predicted by BL97. Unfortunately, the BL97 models never 
exceed densities larger than 3-5 xlO"* cm“^, so no abun¬ 
dance predictions are given for the densities typical of our 
core centers (1-2 xlO^ cm“^. Table 2). This is immaterial 
for CS and CO, since they are so depleted at these densi¬ 
ties that the exact value of the abundance does not affect 
the radiative transfer results. However, it is critical for 
the N 2 H+ and NH 3 analysis, as some BL97 models show 
the beginning of an N 2 H+ and NH 3 abundance decrease 
at the highest densities (most likely due to depletion), and 
this makes any extrapolation to higher densities especially 
unreliable. For those molecules we have explored two pos¬ 
sibilities: that the abundance stays constant at the highest 
densities (unlikely according to BL97 Figs. 2 and 3), and 
that the abundance drops by an order of magnitude for 
densities not studied by BL97 (a conservative depletion 
model). 

Fig. 7 compares the BL97 model predictions for the 
L1517B core with our observations of C^®0, CS, N 2 H+, 
and NH 3 (a similar comparison with LI498 gives the same 
result). As can be seen, the BM94 contraction history 
(dashed lines) predicts C^®0 and CS intensities 5 times 
weaker than observed, while the accretion model predicts 
the correct C^®0 intensity and overestimates the CS emis¬ 
sion by a factor of 3. The weaker C^®0 and CS emission 
in the BM94 model result from its slow core contraction, 
which makes the outer core material rather old (several 
Myr) and, therefore, heavily depleted of these species. The 
accretion model, on the other hand, predicts a younger 
outer core (< 1 Myr) with much less depletion. The higher 
CS emission predicted by the accretion model arises from 
an overestimate of the undepleted CS abundance by at 
least an order of magnitude. 

Both the BM94 and accretion models predict similar 
N 2 H+ and NH 3 abundances, so Fig. 7 only presents the 
results for the accretion model. The higher curves in the 
figure are intensity predictions for the case of no depletion 
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Fig. 7.— Comparison between L1517B observations (solid squares) and expected emission from Bergin & Langer (1997) models (lines). In 
the C^®0 and CS panels, the dashed lines correspond to the BM94 model (their Model 1), and the solid lines indicate the accretion model 
(Model 2). In the N 2 H+ and NH 3 panels, only results for the accretion model are shown. The top lines in these panels represent no depletion, 
and the bottom lines represent a factor of 10 depletion for densities higher than those studied by Bergin & Langer (1997) (see text). All 
models assume dust grains with H 2 O mantles, and to simplify the display, both model and data in the two rightmost panels have been scaled 
down by 0.65 (N 2 H+) and 0.1 (NH 3 ). 


even at the highest densities, while the lower curves as¬ 
sume a moderate depletion of a factor of 10 for densities 
higher than those studied by BL97. As can be seen, the 
models with no depletion agree much better with obser¬ 
vations, although these models are less likely according to 
the BL97 plots. The undepleted N 2 H+ model is partic¬ 
ularly good, and it predicts the correct intensity within 
25% (again, a possible decrease in trapping by hfs was ne¬ 
glected) . The NH 3 prediction, on the other hand, not only 
overestimates the intensity, but it does not reproduce the 
NH 3 abundance increase observed toward the core centers 
(an ortho-para ratio of 2 was assumed). This is a general 
problem of the BL97 models, which predict slight NH 3 
abundance increases at late times, but they are too small 
to fit the NH 3 enhancement inferred from our observations. 
A possible solution is suggested by the recent work of Ne- 
jad & Wagenblast (1999), whose models show a sharp 
enhancement of the NH 3 abundance in the core inner 0.1 
pc while the N 2 H+ fraction remains constant. More de¬ 
tailed modeling of the chemical reactions triggered by the 
disappearance from the gas phase of CO and other species 
is clearly necessary. 

If we take the BL97 models at face value, the previous 
discussion implies that cores condense rapidly and have 
dust grains coated with H 2 O ice. These conclusions, how¬ 
ever, depend critically on the assumed molecular binding 
energies, and they may change as more accurate molecu¬ 
lar parameters are derived. A lower CO binding energy, 
for example, will slow down the rate of CO depletion, 
and this can favor the larger evolution time scales pre¬ 
dicted by ambipolar diffusion. Grain properties, in addi¬ 
tion, may evolve as molecules deplete. An initially H 2 O- 
coated grain may evolve (through CO depletion) into a 
CO-coated grain, and its depletion ability may obviously 
change, or as Caselli et al. (2001c) have proposed, the 
presence of atomic oxygen in the gas phase may maintain 
polar mantles in the densest parts of cores, given that a 
fraction of O will deplete and quickly be transformed in 
H 2 O. These and other uncertainties show that it is prema¬ 
ture to infer core evolution properties from the observed 
pattern of abundances. Still, the reasonable success of the 
BL97 model (and those by Aikawa et al. 2001 and Caselli 
et al. 2001 c) shows that most, or all, the chemical dif¬ 
ferentiation we have observed in starless cores very likely 


arises from depletion and reactions triggered by depletion. 

6.2. A Solution to the NH^/CS Size Discrepancy 

The systematic chemical differentiation of cores we have 
found, and particularly the CS abundance drop, suggests 
a natural explanation for the decade-old discrepancy be¬ 
tween CS and NH 3 observations (e.g.. Fuller & Myers 
1987; Zhou et al. 1989; Pastor et al. 1991; Morata et 
al. 1997). This discrepancy arises because CS rotational 
transitions have larger critical densities than NH 3 inver¬ 
sion transitions. Thus, when mapping a centrally concen¬ 
trated core, one would expect the CS emission to appear 
more peaked and concentrated than the NH 3 emission. In 
practice, however, CS maps are systematically more ex¬ 
tended than NH 3 maps, with ratios between CS and NH 3 
diameters of 1.5 ± 0.3 (Zhou et al. 1989) or « 2 (Myers 
et al. 1991). 

Explanations for the larger extension of the CS emission 
have postulated the scattering of photons by the core outer 
gas (Fuller & Myers 1987), CS abundance enhancement 
by shocks in the outer core layers (Zhou et al. 1989), or 
unresolved clumps (Taylor et al. 1996). However, it has 
not been proved that these mechanisms are truly at work, 
and the NH 3 /CS discrepancy has remained unsolved for 
more than a decade. In this section we explore quantita¬ 
tively how molecular differentiation helps to explain the 
discrepancy, and for this purpose, we use our Monte Carlo 
model results as if they represented observed data and an¬ 
alyze them using standard procedures. 

Our models, by construction, fit simultaneously 
the emission size of C^®O(l-0), CS(2-I), N 2 H+(l- 0 ), 
NH 3 ( 1 , 1 ) and NH 3 ( 2 , 2 ) for each of the 5 cores in our 
sample. This means that these models should naturally 
give rise to any difference in emission size between tracers 
and, therefore, they should reproduce the NH 3 /CS dis¬ 
crepancy. The advantage of reproducing the discrepancy 
with the models is that the internal model parameters are 
well known, so we can use the models as laboratories to 
explore how the discrepancy arises from the core physical 
properties. Taking advantage of our multi-line observa¬ 
tions, we study the behavior of C ^®0 and N 2 H+, in addi¬ 
tion to CS and NH 3 . 

To estimate the sizes of the emission in different 
molecules, we treat our model results as if they were real 
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Fig. 8. — Ratio of map sizes (half-maximum radius) for different molecules as measured from the best-fit Monte Carlo models represented 
in Figs. 3, 5, and 6 . The dashed lines represent the mean value for each molecular pair (see text). Note how the NH 3 sizes are systematically 
smaller than the CS sizes by a factor larger than 2. This systematic size discrepancy, already known but unexplained, is a consequence of the 
central abundance drop of CS, which severely truncates the emission from this molecule and therefore increases the half-maximum radius. 
See text for a complete discussion. 
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data taken with a telescope (FCRAO or 100m Effelsberg), 
and follow the standard procedure of measuring sizes using 
the half maximum radius (Zhou et al. 1989; Myers et al. 
1991). For NH 3 , N 2 H+, and CS this process is easily done 
directly from the model radial profiles, but for C^® 0 , some 
radial profiles present a slight central dip, so the peak in¬ 
tensity does not correspond to the value at the core center. 
For compatibility with the other molecules, we define as 
radius the distance at which the C ^®0 intensity drops to 
half the value at the core center. Had we taken the C^®0 
peak value (off center), the radius would have been larger, 
but only by a very small amount (e.g., 5% in L1517B). 

Figure 8 shows the CS/C^®0, NH 3 /CS, and N 2 H+/CS 
size ratios measured from our Monte Carlo models fol¬ 
lowing the above procedure. As the figure shows, the 
CS radii are consistently smaller than the C^®0 radii 
by about 30%, and the NH 3 and N 2 H+ radii are con¬ 
sistently smaller than the CS radii by a factor of 2 or 
more. Quantitatively, we find that {R{CS)/R{C^^O)) = 
0.73 ± 0.12, {R{NH3)/R{CS)) = 0.39 ± 0.07, and 
{R{N2H+)/R{CS)) = 0.46 ± 0.08. (A direct esti¬ 
mate of these ratios from the sizes of the half max¬ 
imum contours in the maps of Fig. 1 yields 0.77, 
0.39, and 0.54 for R{CS)/R{C^^O), R{NH3) / R{CS) , and 
R{N 2 H'^)/R{CS), respectively, in excellent agreement 
with the values estimated from the Monte Carlo models.) 

The value of {R{CS)/R{C^^O)) in our Monte Carlo 
models is in excellent agreement with the « 0.8 derived by 
Myers et al. (1991) from their sample of cores, and the 
{R{NH 3 )/R{CS)) value also agrees with the « 0.5 mea¬ 
sured by Myers et al. (1991). This good match between 
the relative sizes of the Monte Carlo models and those of 
the general population of cores suggests that the models 
not only reproduce the details of our (small) core sample, 
but that they capture the basic properties of the general 
core population (both starless and with stars). In particu¬ 
lar, the models automatically reproduce the CS/NH 3 size 
discrepancy. 

To understand how the CS/NH 3 size discrepancy arises 
in the models (Fig. 8 central panel), we isolate the factors 
that cause the NH 3 maps to be smaller and the CS maps 
to be larger, and look for the dominant one. The central 


NH 3 abundance enhancement obviously decreases the size 
of the NH 3 maps, and we quantify this effect by comparing 
the radii of the best fit with that of the constant abun¬ 
dance models in Fig. 6 . In this way, we find that best-fit 
models (i.e., abundance enhanced) are about 30% smaller 
than constant-abundance models, which is the right trend 
but is not enough to explain the CS/NH 3 size discrepancy. 
The small effect of the abundance enhancement is consis¬ 
tent with the fact that there is a similar size discrepancy 
between CS and the constant abundance N 2 H+ (see right 
panel in Fig. 8 ). 

Two factors that increase the size of the CS maps are 
optical depth and chemical differentiation. We study the 
effect of optical depth by decreasing the CS abundance 
with a global factor of 22 (to simulate a thin C^^S obser¬ 
vation), and find that this only decreases the average half¬ 
maximum radius by a factor of 20%. This result is consis¬ 
tent with the C^^S observations of L1544 (Fig. 4), which 
show that the rare-isotopomer map is only 25% smaller 
than the CS map, despite the fact that the main isotope 
is extremely thick and self absorbed. We, therefore, con¬ 
clude that optical depth is only a minor factor in the size 
discrepancy. 

Finally, we explore the effect of the central CS abun¬ 
dance drop by comparing the sizes of our best CS models 
with the sizes predicted by constant CS abundance models 
which fit the observed central intensity. These constant- 
abundance models are a full factor of 2 smaller than the 
models with abundance drop, so they have sizes more com¬ 
parable to those of NH 3 and N 2 H+. This means that the 
major contributor to the CS/NH 3 model discrepancy is 
the central drop in CS abundance, which decreases the 
peak CS intensity by a factor of several (see Fig. 4 bot¬ 
tom panel), and this in turn increases the radius of the 
half-maximum contour by about a factor of 2 . 

6.3. The NH 3 /CS Linewidth Discrepancy and the 
Linewidth-Size Relation in Starless Cores 

A problem commonly associated with the NH 3 /CS size 
discrepancy is the systematic difference in linewidths be¬ 
tween these two dense gas tracers. In their comparison 
between NH 3 and CS observations of dense cores, Zhou et 
al. (1989) found that AV{CS)/AV{NH 3 ) = 2.0 ± 0.6, 
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Fig. 9.— Comparison of CS and NH 3 linewidths from dense cores. Crosses are data from Zhou et al. (1989) (starless cores only), and 
filled squares are results from our Monte Carlo models. Both real data and models lie below the line of equal CS and NH 3 linewidth, implying 
that CS lines are systematically broader than NH 3 lines. Although this linewidth discrepancy has been traditionally interpreted as resulting 
from the systematic increase of gas turbulence with core radius, none of our models has this property. Optical depth broadening and self 
absorption are the two causes of the broader CS lines in the Monte Carlo models (see text). 


which is too large to be explained by optical-depth broad¬ 
ening of the CS lines. Given the additional larger size of 
the CS emission, the traditional interpretetaion of this ef¬ 
fect has been in terms of a systematic increase in the gas 
turbulence with radius, following the so-called linewidth- 
size relation (e.g.. Fuller & Myers 1992). 

Given the above explanation, it appears surprising that 
all our cores (but L1544) are well fit by models with con¬ 
stant turbulent velocity, and that for L1544, our fit re¬ 
quires a decrease of linewidth with radius, not an increase. 
This result is significant because the line modeling is very 
sensitive to any linewidth variation along the line of sight 
through the simultaneous fit of the CS(2-1) self absorption 
(which arises from gas at large radius) and the fit of the 
N 2 H+ and NH 3 lines (which represent gas from the inner¬ 
most core). The success of the constant-linewidth models 
suggests that, in most cores, there are no systematic vari¬ 
ations in the turbulent-velocity field from the core center 
up to radii of at least 0.1 pc (core edge). 

To study whether our modeling is consistent with the 
NH 3 /CS linewidth discrepancy observed in other cores, 
we again use our model data as if they were real obser¬ 
vations, and measure CS and NH 3 linewidths from the 
predicted spectra. To do this, we first convert the model 
spectra into CLASS format and then use that program 
to fit the lines. Following standard procedure, we fit the 
CS(2-1) lines with simple Gaussians, and the NHs)!,!) 
multiplets with the standard hyperfine analysis. We apply 
this procedure even to the obviously self-absorbed L1544 
CS spectrum, and in Fig. 9 plot our model data (black 
squares) together with the data from Zhou et al. (1989) 
for starless cores (crosses). No attempt has been made to 
correct our model data from instrumental broadening, as 
the resolution of our model ensures that this effect is less 
than 5%. 

As Fig. 9 shows, our model data lie in the same re¬ 


gion where we find the data from observations of other 
cores, i.e. consistently below the line of equal NH 3 and 
CS linewidth. This shows that our models naturally dis¬ 
play the NH 3 /CS linewidth discrepancy, the way they did 
for the size discrepancy. The model data, in addition, 
predict AV{CS)/AV{NH 3 ) = 1.5 ± 0.4, which is consis¬ 
tent with the ratio derived by Zhou et al. (1989). The 
clustering of most of our points in the region of smallest 
linewidth is probably due to the fact that we selected the 
most quiescent cores possible (Section 2), and in fact, all 
cores but L1544 were modeled using the same turbulent 
speed (L1544 is the outlying point with large AV{CSj). 
From Fig. 9 and the measured linewidth ratio, we con¬ 
clude that a significant NH 3 /CS linewidth discrepancy can 
be produced even in the presence of a constant turbulence 
velocity field, and that it does not require the presence of 
a systematic increase of linewidth with size. 

To understand how this NH 3 /CS linewidth discrepancy 
arises in constant turbulence gas, we have analyzed how 
lines are generated in our models and have found two main 
factors that contribute to the apparent increase of the CS 
linewidth over that of NH 3 . First, there is a bias in the 
way the two linewidths are measured: in the standard 
procedure, CS linewidths are derived by fitting a single 
Gaussian while NH 3 linewidths are estimated with a hy¬ 
perfine fit. This method has the somewhat perverse effect 
of correcting he usually not-very-thick NH 3 line for opacity 
broadening, while not correcting at all the always-thick CS 
spectra. As Zhou et al. (1989) mention, this effect con¬ 
tributes to, but cannot be the only factor in, the NH 3 /CS 
linewidth discrepancy. 

The additional factor that enhances the CS linewidth 
is the presence of self absorption. Self absorption par¬ 
tially truncates the CS line, decreasing its peak intensity 
and systematically increasing the size of its half-maximum 
width. This effect is similar to the way the depletion 
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of CS at the core center truncates the map of emission 
and increases its half-maximum radius. Together with the 
linewidth-measurement bias, it is responsible for the con¬ 
sistently larger CS linewidths in our models (systematic 
velocity fields have only a minor effect). 

The effect of self absorption in the final linewidth is 
very sensitive to the velocity displacement between the 
emitting and absorbing gas. A large-enough displacement 
(like in L1498 and L1517B, see Fig. 3) produces a single- 
peaked line slightly broader than the NH 3 line, while no 
displacement (like in L1544) produces a double-peaked line 
which yields an extremely large linewidth when fitted with 
a Gaussian. This effect shifts the L1544 point very far to 
the right in Fig. 9, and it suggests that the points with 
large AC(C'S') from Zhou et al. (1989) represent spectra 
with deep central self-absorptions, although further work 
is needed to clarify this issue. (Note that the strongly 
double-peaked CS line from L1544 would have appeared 
rather diluted in the 0.09 km s“^ channels of Zhou et al. 
1989, and that these authors recognize that some of their 
lines are truly self-absorbed.) 

Our conclusion that the gas turbulence does not increase 
with radius in the inner part of our cores, based on the 
analysis of line-of-sight motions, agrees with recent stud¬ 
ies of turbulent motions in the plane of the sky for other 
cores. Using NH 3 as a single tracer, Barranco & Goodman 
(1998) find negligible linewidth variations up to distances 
of 0.1 pc, which is inconsistent with the factor-of-3 in¬ 
crease expected from the standard linewidth-size relation 
between 0.01 and 0.1 pc. This suggests that the common 
interpretation of a turbulence increase in the inner 0.1 pc 
of some cores, based on the combination of NH 3 and GS 
data, may arise from the combination of optical depth ef¬ 
fects (which increase the CS linewidth) and a central abun¬ 
dance drop (which increases the CS half-maximum radius). 

What happens at radii larger than 0.1 pc is still unclear. 
Large-scale studies show that there is a true linewidth- 
size relation when considering whole clouds (e.g., Larson 
1981), and Goodman et al. (1998) have suggested that 
turbulence starts to increase at radii larger than 0.1 pc. 
Unfortunately their NH 3 data are not sensitive enough to 
follow this trend very far. An alternative view suggested 
by some of our CS and C^®0 data is that the large-scale 
broadening may result from the overlap of different cloud 
fragments moving at different velocities. Further work on 
large-scale modeling of molecular lines should be carried 
out to clarify this crucial topic. 

6.4. Further Consequences of Chemical Differentiation 

The chemical differentiation of starless cores has addi¬ 
tional implications for studies of low-mass star formation, 
and in this section we briefly discuss some aspects of this 
problem. CS, for example, is one of the molecules of choice 
for density determinations (e.g., Evans 1999), thus, its 
vanishing from the gas phase at densities of a few 10 "^ cm“^ 
implies that any attempt to measure with this molecule 
the density profile of a concentrated core seems doomed 
to miss the inner core region and, therefore, to underes¬ 
timate the central density. This same problem may af¬ 
fect other classical density tracers like HCO'*', HCN, and 
H 2 CO (see BL97 for model predictions), and a system¬ 
atic analysis is in progress for these and other molecules 


(Tafalla et ah, in preparation). The constant abundance 
of N 2 H+ suggests that this molecule may be a robust den¬ 
sity tracer (see Caselli et al. 1999), but its complicated 
hyperfine structure and a lack of accurate collisional rate 
coefficients make its radiative transfer much more difficult 
to solve. Although it may be argued that dust contin¬ 
uum observations may suffice to derive core density profiles 
(our approach in section 4.2), confirmation of these esti¬ 
mates with molecular observations seems necessary, given 
the present uncertainties in the dust emission properties. 

An additional issue in the study of starless cores poten¬ 
tially affected by chemical differentiation is the search for 
infall motions. Infall searches are usually done combining 
observations of optically thick and thin tracers, often from 
different molecular species (e.g., Lee et al. 1999). A com¬ 
mon choice consists of GS (thick) and N 2 H+ (thin), two 
species which we have shown to reside in different parts 
of the core. This disjoint distribution of tracers compli¬ 
cates the interpretation of the spectral signatures and can 
only be avoided using isotopomer pairs of different opti¬ 
cal depth. The sharp central abundance drop of CS, in 
addition, implies that infall searches using this molecule 
are only sensitive to the motions of the core’s outer lay¬ 
ers. Thus, to probe motions close to the core center, 
a molecular species resistant to differentiation/depletion, 
like N 2 H“'', is needed. The lower optical depth of the N 2 H“'' 
lines, however, complicates detecting inward motions with 
the standard self-absorption analysis (but see Williams et 
al. 1999 for a claim of N 2 H+ self-absorption in L1544). 
The search for high velocity wings in the lines of this tracer 
appears to be a complementary and promising technique 
(Bourke et al. 2001). 

Chemical differentiation does not only represent an 
added difficulty in the study of starless cores. The time 
dependence of the underlying chemical processes suggests 
that detailed chemical modeling should allow us to recon¬ 
struct core contraction histories from observed depletion 
patterns. We cautioned in section 6.1 that the fast con¬ 
traction scenario derived from the BL97 model is still ten¬ 
tative, as it depends on the exact value of the molecu¬ 
lar binding energies, and that an accurate determination 
of these critical parameters is still needed to definitively 
test the longer time scales predicted by ambipolar diffu¬ 
sion models. Before such quantitative time estimates are 
possible, we can compare the abundance patterns in dif¬ 
ferent cores in an attempt to construct a qualitative time 
sequence of core evolution. Most cores in our sample have 
similar central densities (~ 10 ® cm“®) and similar deple¬ 
tion patterns (Tables 4 and 5), so it seems likely that they 
are at a similar evolutionaly stage (which may be just an 
effect of our selection criteria). Interestingly enough, our 
sample contains cores with both inward and outward mo¬ 
tions (LI498 and L1517B, respectively), suggesting that 
the presence of these large-scale velocity patterns has lit¬ 
tle to do with the core evolutionary state. 

The only clue to an evolutionary pattern in our sample 
comes from the most extreme object, LI544. This core has 
density higher than the others by an order of magnitude 
(~ 10 ® cm“®), and its abundance pattern slightly differs 
in the sense that the abundance variations of CS, C®®0, 
and NH 3 occur at higher densities than in the other cores 
(Tables 4 and 5). This could mean that a very rapid con- 
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traction has kept frozen a chemical pattern characteristic 
of lower densities or that the core has contracted faster 
than the others, which would not allow enough time for 
the moderately dense gas to be chemically processed. 

Clearly additional work is needed before going further 
with the above speculations. The chemical composition 
of cores with different densities and conditions should be 
measured in detail. Further modeling of the chemical evo¬ 
lution of cores as they contract is also needed. Finally, it 
is important to identify key molecules that allow us to dis- 
tiguish between different chemical scenarios, and to char¬ 
acterize the chemical behavior of standard high-density 
tracers. The potential benefits of understanding the chem¬ 
ical evolution of cores make exploring these and related 
questions an urgent enterprise. 

7. CONCLUSIONS 

We have mapped 5 starless cores (L1498, L1495, 
L1400K, L1517B, and L1544) in Ci®O(l-0), CS(2-1), 
N 2 H+( 1 - 0 ), NH 3 ( 1 , 1 ) and (2,2), Ci^O(l-O) (3 cores), 
C^‘*S(2-1) (1 core), and the 1.2 mm continuum (4 cores), 
and complemented these data with a published 1 . 2 mm 
continuum map of L1544 (Ward-Thompson et al. 1999). 
For each core, we have self-consistently determined the ra¬ 
dial profile of density (from the 1 . 2 mm continuum), tem¬ 
perature (from NH 3 ), and molecular abundance (using a 
radiative transfer Monte Carlo model). To do this, we 
have fit simultaneously the central spectrum for each line 
and the radial profile of integrated intensity. As a result 
of this work, we conclude the following: 

1. Each core is well fit with a radial density profile of the 
form n{r) = no/(l -f (?"/ro)“) where no ranges from about 
10® to 10® cm“®, ro is of the order of 3000-10000 AU, and 
a ranges from 2 to 4. This type of profile naturally fits 
the central flattening and the large-r power-law behavior 
found by previous authors. 

2. All cores show evidence for a strong central drop 
in the abundance of C^®0 and CS. This drop can be fit¬ 
ted with a negative exponential dependence on density 
(exp(—n(r)/n(i)), with an e-folding Ud parameter of 2-6 
X lO'^ cm“®. C^^O(l-O) and C®^S(2-1) observations of sev¬ 
eral cores confirm these results and show that our radiative 
transfer models can naturally fit this rare isotopomer emis¬ 
sion without extra free parameters. Although limited by 
the resolution of our data (~ 50"), we estimate that the 
C^®0 and CS central abundance drops are at least of one 
or two orders of magnitude in each of the cores. 

3. In contrast to C^®0 and CS, the N 2 H+ abundance 
seems to be constant inside each core. The combination 
of narrow N 2 H+(l- 0 ) central spectra and the presence of 
narrow CS(2-1) self absorptions indicate that the turbu¬ 
lent linewidth in each core is constant with radius. 

4. The NH 3 ( 1 , 1 ) and (2,2) lines are well fit with con¬ 
stant temperature models of 9.5 K (8.75 K in L1544) in 
which the NH 3 abundance increases toward the core cen¬ 
ter. The fact that the temperature is constant at densities 
for which gas and dust may be thermally coupled supports 


the assumption of constant dust temperature used in our 
1 . 2 mm continuum analysis. 

5. The combination of a central C^®0 and CS abun¬ 
dance drop, a constant N 2 H''' abundance, and a central 
NH 3 abundance increase seems to be a systematic charac¬ 
teristic of starless cores. This suggests that star-forming 
material is very chemically inhomogeneous before starting 
to collapse, and that searches for infall in starless cores 
should take this fact into account. 

6 . The pattern of chemical differentiation, especially for 
C^®0, CS, and N 2 H+ is in qualitative agreement with de¬ 
pletion models. Quantitative agreement between the data 
and these models can be achieved for certain choices of 
their parameters, which suggest the presence of polar ices 
on grains and a history of core contraction more rapid 
than predicted by subcritical ambipolar-diffusion models. 
These models, however, seem to overestimate N 2 H+ de¬ 
pletion at high densities and to miss the observed central 
NH 3 abundance increase. 

7. The presence of chemical differentiation in cores auto¬ 
matically explains the systematic difference between map 
sizes of different tracers, a problem which has remained un¬ 
explained for more than a decade. Detailed Monte Carlo 
radiative transfer models show that CS maps, because of 
the central abundance decrease and radiative transfer ef¬ 
fects, are expected to appear at least twice as large as 
NH 3 and N 2 H+ maps. These latter molecules are there¬ 
fore more faithful tracers of the dense core material. 

8 . Comparing the CS and NH 3 linewidths from our 
Monte Carlo radiative transfer calculations, we find a sys¬ 
tematic NH 3 /CS linewidth discrepancy similar to that 
found by other authors (e.g., Zhou et al. 1989). These 
linewidth discrepancies have been traditionally interpreted 
as resulting from a systematic increase in the turbulent 
linewidth with radius, but our models show that it can 
arise from a combination of optical depth and self absorp¬ 
tion in the CS line. This suggests that the linewidth-size 
relation does not apply in the inner 0.1 pc of cores. 
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Table 1 


Sample Cores. 


Core 

a(1950) 

h m s 

(5(1950) 

O / // 

L1498 

04 07 50.0 

25 02 13 

L1495 

04 11 02.7 

28 00 43 

L1400K 

04 26 51.0 

54 45 27 

L1517B 

04 52 07.2 

30 33 18 

L1544 

05 01 14.0 

25 07 00 


Table 2 

Density fits from continuum data*-®^ 


Core 

Center*^^^ 

("/') 

T^ic) 

(K) 

no/lO^ 
(°) (cm“^) 

ro 

n 

a 

L1498 

(0,0) 

9.5 

0.6 

-40 

1.0 

75 

4 

L1495 

(20,40) 

9.5 

1.0 

0 

1.1 

45 

2 

L1400K 

(-20,-40) 

9.5 

0.3 

45 

1.2 

30 

2 

L1517B 

(-10,-20) 

9.5 

1.0 

0 

2.2 

35 

2.5 

L1544 

(-13,-21) 

8.75 

0.6 

-45 

14 

20 

2.5 


= no/[l + (r/ro)“] 

Offset with respect to coordinates in Table 1. 

^‘^^Dust temperature, assumed equal to gas temperature (see 
text for details). 

Aspect ratio for elliptical average. 

Position angle for elliptical average. 


Table 3 


Core velocity parameters. 


Core 

Vlsr 
( km s“^) 

AV 

(km s“^) 

dV/dr 

(km s“^ pc“^) 

L1498 

7.80 

0.17(“) 

-0.47 

L1495 

6.80 

0.17(“) 

0.93 

L1400K 

3.16 

0.17(“) 

1.55 

L1517B 

5.78 

0.17(“) 

0.78 

LI 544 

7.20 

0.20^ 

0.00 


Spatially constant. 
Variable, see section 5.1. 
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Table 4 

Best-Fit Model Parameters for and CS. 


CORE 

^0 

Ci ®0 

Ud (cm“^) 

fib) 

^0 

CS 

Ud (cm-3) 

fib) 

L1498 

1.7 10"^ 

2.0 lO"* 

2 10-2 

3 10-9 

3.0 104 

8 10-2 

L1495 

1.7 10-^ 

2.2 lO"* 

4 10-2 

5 10-9 

1.5 104 

8 10-3 

L1400K 

1.7 10-^ 

1.4 10^ 

5 10-3 

9 10-9 

1.5 104 

7 10-3 

L1517B 

1.7 10-^ 

2.1 10 ^ 

4 10-4 

3 10-9 

4.0 104 

2 10-2 

L1544 

1.7 10-^ 

5.5 10^ 

1 10-3 

3 10-9 

1.7 103 

2 10-2 


(“)X(r) = Xoexp(-n(r)/nd) 

Model abundance drop: X{r = 20”)/X{r = 100") 


Table 5 

Best-Fit Model Parameters for N2H+ and NH3. 


CORE 

N 2 H+ 

Xo(3) 

NH 3 

/3(3) 

fic) 

L1498 

9.0 10 -“ 

1.0 10-3 

1.0 

4.1 

L1495 

2.4 10-4° 

2.3 10-3 

1.0 

4.8 

L1400K 

1.3 10-49 

4.0 10-9 

0.0 

1.0 

L1517B 

9.5 10-44 

1.7 10-3 

1.0 

12 

LI 544 

7.5 10-44 

4.0 10-9 

0.3 

2.7 


Constant N 2 H“*' abundance. 

('')X(NH3) = Xoin/no)^ 

Model abundance enhancement: X{r = 
20")/X{r = 100"). 



